library("knitr") # for knitting stuff
library("kableExtra") # for markdown tables
library("lme4") # for linear mixed effects models
library("broom.mixed") # for tidying mixed models results
library("brms") # Bayesian regression with Stan
library("corrr") # for tidy correlation matrix
library("xtable") # for latex tables
library("jsonlite") # for reading json files
library("janitor") # cleans stuff
library("Hmisc") # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png") # adding pngs to images
library("grid") # functions for dealing with images
library("plotly") # 3D scatter plot
library("egg") # for geom_custom
library("clValid") # for validating clustering solutions
# library("ggtern") # for ternary plots
library("patchwork") # for figure panels
library("tidyverse") # for wrangling, plotting, etc. # setting some knitr options
opts_chunk$set(comment = "",
fig.show = "hold")
# set the ggplot theme
theme_set(theme_classic())
# suppress summarise() grouping warning
options(dplyr.summarise.inform = F)
# function for printing out html or latex tables
print_table = function(data, format = "html", digits = 2){
if(format == "html"){
data %>%
kable(digits = digits) %>%
kable_styling()
}else if(format == "latex"){
data %>%
xtable(digits = digits) %>%
print(include.rownames = F,
booktabs = T)
}
}
# root mean squared error
rmse = function(x, y){
return(sqrt(mean((x - y)^2)))
}
actual_counterfactual_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = clipindex,
y = value,
fill = colorindex,
group = model)) +
stat_summary(geom = "bar",
fun = "mean",
color = "black",
position = position_dodge(0.7),
width = 0.7) +
stat_summary(geom = "linerange",
fun.data = "mean_cl_boot",
position = position_dodge(0.7),
size = 1) +
geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
shape = 21,
color = "black",
size = 1,
alpha = 0.2) +
scale_fill_manual(values = c("red2", "green2", "black")) +
facet_grid(index_actual ~ index_counterfactual) +
geom_text(data = df.text %>%
drop_na(label),
mapping = aes(x = x, y = y, label = as.character(label)),
size = 6,
position = position_dodge(width = .7),
hjust = 0.5) +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 105),
expand = F,
clip = "off") +
labs(y = ylabel) +
theme(text = element_text(size = 20),
legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
panel.spacing.y = unit(0.8, "cm"),
plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
axis.title.x = element_blank(),
panel.grid = element_blank(),
strip.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
# print(p)
}
actual_counterfactual_threeballs_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = ball,
y = value,
group = model,
fill = colorindex)) +
stat_summary(fun = "mean",
geom = "bar",
color = "black",
position = position_dodge(0.9),
width = 0.9) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
color = "black",
position = position_dodge(0.9)) +
geom_point(data = x %>%
mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(jitter.width = 0.3,
jitter.height = 0,
dodge.width = 0.9),
shape = 21,
color = "black",
size = 1,
alpha = 0.15) +
facet_wrap(~clip, ncol = 8) +
geom_text(data = df.text,
mapping = aes(x = x,
y = y,
label = as.character(label)),
size = 5,
position = position_dodge(width = .9),
hjust = 0.5,
fontface = "bold") +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 120),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25)) +
scale_fill_manual(values = c("red2", "green2", "black")) +
scale_color_manual(values = c("red2", "green2", "black")) +
labs(y = ylabel) +
theme(text = element_text(size = 16),
legend.position = "none",
axis.title.x = element_blank(),
panel.spacing.y = unit(0.3, "cm"),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
}# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv",
header = T,
stringsAsFactors = F)
# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp1.clipinfo = tibble(clip = 1:18,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 1),
outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1,
0, 1, 1, 0, 0, 1, 0, 1, 1),
index_actual = rep(c("actual miss",
"actual close",
"actual hit"), each = 6),
index_counterfactual = rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3)) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp1.causal.long = df.exp1.causal$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")
df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")# counterfactual condition
df.exp1.counterfactual %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 32 | 8.8 | 19 | 7.5 | 3.12 |
# causal condition
df.exp1.causal %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 33 | 10.3 | 21 | 9.65 | 6.58 |
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]
df.exp1.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp1.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup() %>%
left_join(df.exp1.clipinfo, by = "clip")
# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp1.counterfactual.model = df.exp1.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp1.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)df.exp1.causal.model = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(empirical = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.counterfactual.model %>%
select(-c(noise, sse)),
by = "clip") %>%
rename(simulation = prediction) %>%
mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))set.seed(1)
df.plot = df.exp1.counterfactual.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.)/2)) %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model, levels = c("rating", "prediction"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
# width = 8,
# height = 6)df.exp1.counterfactual.means %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 1 | 0.97 | 10.95 |
set.seed(1)
model.name = c("empirical")
# model.name = c("simulation")
df.plot = df.exp1.causal.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, simulation, empirical),
names_to = "model",
values_to = "value") %>%
filter(model %in% c(model.name, "rating")) %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(
colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index.actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index.counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)
# ggsave("../../figures/plots/exp1_causal_bars.pdf",
# width = 8,
# height = 6)df.exp1.causal.long %>%
mutate(rating = abs(rating)) %>%
group_by(clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(empirical_r = cor(rating, empirical),
empirical_rmse = rmse(rating, empirical),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.96 | 8.57 | 0.93 | 15.15 |
df.exp1.causal.posterior = df.exp1.causal.long %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp1_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)
func_posterior_difference = function(data, clip1, clip2){
data %>%
mutate(difference = .data[[clip1]] - .data[[clip2]]) %>%
pull(difference) %>%
mean_hdci() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>%
print()
}
func_posterior_difference(data = df.exp1.causal.posterior,
clip1 = "8",
clip2 = "13") difference
1 (0.26 [-6.51, 7.31])
fit.brm_exp1_causal %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | -18.33 | 3.36 | -24.86 | -11.66 |
| fixed | cond | index_actualactualclose | 4.31 | 3.40 | -2.54 | 10.94 |
| fixed | cond | index_actualactualhit | 4.06 | 4.84 | -5.51 | 13.27 |
| fixed | cond | index_counterfactualcounterfactualclose | -17.52 | 3.37 | -24.36 | -10.84 |
| fixed | cond | index_counterfactualcounterfactualhit | -61.37 | 4.38 | -69.95 | -52.89 |
| fixed | cond | outcome_actual | 92.36 | 5.11 | 82.21 | 102.38 |
df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
header = T,
stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp2.clipinfo = tibble(clip = 1:20,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
index_actual = c(rep(c("actual miss",
"actual close",
"actual hit"),
each = 6),
rep("practice", 2)),
index_counterfactual = c(rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3),
rep("practice", 2))) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp2.causal_first.long = df.exp2.causal_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("causal", "counterfactual"), each = 20),
max(participant)),
condition = "causal_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("counterfactual", "causal"), each = 20),
max(participant)),
condition = "counterfactual_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
bind_rows(df.exp2.counterfactual_first.long %>%
mutate(participant = participant +
max(df.exp2.causal_first.long$participant)))# counterfactual condition
df.exp2.counterfactual_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time) %>%
bind_rows(df.exp2.causal_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time)) %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 82 | 35 | 12.2 | 45 | 21.25 | 5.11 |
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]
df.exp2.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp2.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp2.counterfactual.means = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
# mutate(rating = abs(rating)) %>%
group_by(clip) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.clipinfo, by = "clip")
# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp2.counterfactual.model = df.exp2.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
group_by(clip, condition) %>%
summarize(rating = mean(rating)) %>%
pivot_wider(names_from = condition,
values_from = rating) %>%
left_join(df.exp2.combined.long %>%
filter(question == "counterfactual", clip <= 18) %>%
group_by(clip) %>%
summarize(combined = mean(rating)),
by = "clip") %>%
left_join(df.exp2.counterfactual.model,
by = "clip") %>%
select(-c(sse, noise)) %>%
rename(simulation = prediction) %>%
mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
.funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))set.seed(1)
df.plot = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model,
levels = c("rating", "prediction")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
# width = 8,
# height = 6)df.exp2.counterfactual.means %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating, prediction),
simulation_rmse = rmse(rating, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 0.9 | 0.96 | 13.46 |
# random intercepts & slopes model
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual) +
(1 + index_actual +
index_counterfactual | participant),
data = df.exp2.combined.long %>%
na.omit() %>%
filter(question == "counterfactual"),
cores = 2,
seed = 1,
file = "cache/brm_exp2_counterfactual")fit.brm_exp2_counterfactual %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | 13.38 | 2.58 | 8.30 | 18.38 |
| fixed | cond | conditioncounterfactual_first | -1.11 | 3.65 | -7.99 | 6.04 |
| fixed | cond | index_actualactualclose | -1.48 | 2.78 | -6.79 | 3.95 |
| fixed | cond | index_actualactualhit | -4.42 | 2.40 | -9.15 | 0.20 |
| fixed | cond | index_counterfactualcounterfactualclose | 37.94 | 3.33 | 31.43 | 44.60 |
| fixed | cond | index_counterfactualcounterfactualhit | 73.80 | 4.17 | 65.54 | 81.99 |
| fixed | cond | conditioncounterfactual_first:index_actualactualclose | -0.55 | 3.84 | -8.29 | 6.79 |
| fixed | cond | conditioncounterfactual_first:index_actualactualhit | -0.85 | 3.41 | -7.60 | 5.87 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualclose | 1.23 | 4.79 | -8.37 | 10.75 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualhit | -0.02 | 5.92 | -11.96 | 11.59 |
set.seed(1)
func_exp2_causal_plot = function(condition.name, model.name){
df.plot = df.exp2.combined.long %>%
filter(question == "causal", clip <= 18) %>%
filter(condition %in% condition.name) %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c("rating", all_of(model.name)),
names_to = "model",
values_to = "value") %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
actual_counterfactual_plot(df.plot, "causal judgment")
}
plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
model.name = "combined"),
func_exp2_causal_plot(condition.name = "causal_first",
model.name = "combined"))
wrap_plots(plot.list, ncol = 2)
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# model.name = "combined"
#
# func_exp2_causal_plot(condition.name = condition.name,
# model.name = model.name)
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
# width = 8,
# height = 6)Figure 3.1: Left panel: Counterfactual judgments first; right panel: Causal judgments first
df.data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
mutate(rating = abs(rating)) %>%
group_by(condition,
clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual"))
df.data %>%
group_by(condition) %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined)) %>%
print_table()| condition | empirical_r | empirical_rmse |
|---|---|---|
| causal_first | 0.92 | 13.78 |
| counterfactual_first | 0.99 | 5.27 |
df.data %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.95 | 10.43 | 0.92 | 19.09 |
# random intercepts + slopes model
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual + outcome_actual) +
(1 + index_actual + index_counterfactual +
outcome_actual | participant),
data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18),
cores = 2,
seed = 1,
file = "cache/brm_exp2_causal")df.exp2.causal.posterior = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
group_by(clip, index_actual, index_counterfactual,
outcome_actual, condition) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp2_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, condition, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)
# difference between clips 13 and 14 versus 8
df.exp2.causal.posterior %>%
mutate(difference = (`13` + `14`) / 2 - `8`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-0.55 [-12.09, 10.97])
# difference between clips 5 and 6 versus 11
df.exp2.causal.posterior %>%
mutate(difference = (`5` + `6`) / 2 - `11`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-5.48 [-13.8, 2.73])
fit.brm_exp2_causal %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
# print_table(format = "latex")
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | -24.07 | 3.50 | -30.87 | -17.24 |
| fixed | cond | conditioncounterfactual_first | 11.47 | 4.91 | 1.62 | 21.22 |
| fixed | cond | index_actualactualclose | 8.17 | 3.44 | 1.47 | 14.79 |
| fixed | cond | index_actualactualhit | 13.42 | 4.90 | 3.74 | 22.69 |
| fixed | cond | index_counterfactualcounterfactualclose | -30.20 | 3.72 | -37.61 | -23.19 |
| fixed | cond | index_counterfactualcounterfactualhit | -50.29 | 4.60 | -59.57 | -41.47 |
| fixed | cond | outcome_actual | 98.53 | 5.40 | 87.81 | 108.98 |
| fixed | cond | conditioncounterfactual_first:index_actualactualclose | -5.39 | 4.94 | -15.00 | 4.35 |
| fixed | cond | conditioncounterfactual_first:index_actualactualhit | -16.97 | 7.02 | -30.40 | -3.33 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualclose | -8.13 | 5.17 | -18.18 | 2.17 |
| fixed | cond | conditioncounterfactual_first:index_counterfactualcounterfactualhit | -25.51 | 6.44 | -38.16 | -12.98 |
| fixed | cond | conditioncounterfactual_first:outcome_actual | 10.72 | 7.60 | -4.16 | 26.02 |
# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
clipindex = rep(1:2, 16))
# COUNTERFACTUAL JUDGMENTS
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty",
"smoothness", "time", "condition"),
n()/6),
participant = rep(1:(n()/6), each = 6)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, condition, gender, age, time, smoothness, difficulty)
# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "rating"), n()/2),
participant = rep(1:(n()/64), each = 64),
order = rep(rep(1:32, each = 2), n()/64)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.counterfactual.demographics %>%
select(participant, ball = condition),
by = "participant") %>%
select(participant, ball, clip, everything()) %>%
arrange(participant, clip)
# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.causal.demographics = df.exp3.causal$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty", "smoothness", "time",
"counterbalance", "replayType", "experiment"),
n()/8),
participant = rep(1:(n()/8), each = 8)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, gender, age, counterbalance, time, smoothness, difficulty)
# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
participant = rep(1:(n()/(32*6)), each = (32*6)),
order = rep(rep(1:32, each = 6), n()/(32*6))) %>%
filter(!name %in% c("x", "y", "z")) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.causal.demographics %>%
select(participant, counterbalance),
by = "participant") %>%
mutate(A = ifelse(counterbalance == 1, rating1, rating2),
B = ifelse(counterbalance == 1, rating2, rating1)) %>%
select(-rating1, -rating2) %>%
pivot_longer(cols = c(A, B),
names_to = "ball",
values_to = "rating") %>%
select(participant, clip, ball, rating, order) %>%
arrange(participant, clip, ball)# counterfactual condition
df.exp3.counterfactual.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 80 | 33 | 10.1 | 34 | 18.08 | 4.63 |
# causal condition
df.exp3.causal.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 34 | 10.5 | 21 | 21.19 | 4.96 |
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]
df.exp3.model = files %>%
map_dfr(~ read_csv(str_c(path, .))) %>%
rename(clip = trial)# calculate mean counterfactual judgments
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
# find noisy simulation model that best predicts the mean counterfactual
# judgments by minimizing the sum of squared errors
df.exp3.model.counterfactual = df.exp3.model %>%
group_by(noise) %>%
select(clip, contains("whether"), noise) %>%
pivot_longer(cols = c(A_whether, B_whether),
names_to = "ball",
values_to = "prediction") %>%
mutate(ball = str_remove(ball, "_whether")) %>%
arrange(noise, clip, ball) %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp3.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)# calculate mean causal judgments
df.exp3.causal.means = df.exp3.causal.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
df.exp3.model.causal = df.exp3.model %>%
# take into account difference-making
mutate_at(.vars = vars(A_how:A_robust),
.funs = list(~ . * A_difference)) %>%
mutate_at(.vars = vars(B_how:B_robust),
.funs = list(~ . * B_difference)) %>%
# choose model based on best fit with counterfactual data
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(clip, ball:robust)l.features = fromJSON("data/features.json")
df.features.initial = l.features[["appearance"]] %>%
cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
"ballA_velx",
"ballA_vely",
"ballB_velx",
"ballB_vely",
"ballE_velx",
"ballE_vely")) %>%
mutate(clip = 1:nrow(.)) %>%
pivot_longer(cols = starts_with("ball")) %>%
separate(name, into = c("ball", "property")) %>%
pivot_wider(names_from = property, values_from = value) %>%
mutate(ball = str_remove_all(ball, pattern = "ball"))
# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")
for (i in 1:32) {
df.tmp = data.frame()
for (j in 1:length(ballnames)) {
ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
if (ncollisions > 0) {
tmp = data.frame(
ball = ballnames[j],
object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
ncollision = 1:ncollisions
)
df.tmp = df.tmp %>%
rbind(tmp)
}
}
df.features.collisions = df.features.collisions %>%
rbind(df.tmp %>%
mutate(clip = i) %>%
select(clip, ball, everything()))
}
# find end of clip
tmp = df.features.collisions %>%
group_by(clip) %>%
filter(ball == "ballE",
str_detect(object, "goal|Left|Right")) %>%
group_by(clip) %>%
mutate(endclip = max(time)) %>%
select(clip, endclip) %>%
ungroup() %>%
distinct()
# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
left_join(tmp, by = "clip") %>%
mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
group_by(clip) %>%
filter(time <= endclip)
# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
filter(ball == "E") %>%
mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
select(clip, E.moving)
# contact with ball E
tmp.contactE = df.features.collisions %>%
filter(object == "ballE") %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
mutate(contactE = 1) %>%
select(clip, ball, contactE) %>%
group_by(clip, ball) %>%
summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate(contactE = ifelse(is.na(contactE), 0, contactE))
# number of collisions
tmp.ncollisions = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
group_by(clip, ball) %>%
count() %>%
rename(ncollision = n)
# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
count(clip, ball, object) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
group_by(clip) %>%
summarize(AE = any(AE %>% as.logical()) * 1,
BE = any(BE %>% as.logical()) * 1) %>%
mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
select(clip, A, B) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "exclusive")
# impact
func_angle = function(x, y) {
dot.prod = x %*% y
norm.x = norm(x, type = "2")
norm.y = norm(y, type = "2")
theta = acos(dot.prod / (norm.x * norm.y))
as.numeric(theta)
}
# impact on ball E
tmp.impactE = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ball"),
ball == "ballE") %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely),
c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely,
post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(E.speed.diff = speed.diff,
E.direction.diff = direction.diff)
# total impact
tmp.impactTotal = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ballA|ballB")) %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely),
c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely,
post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball,
levels = c("ballA", "ballB"),
labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(total.speed.diff = speed.diff,
total.direction.diff = direction.diff)
# transfer of force
tmp.transfer = df.features.collisions %>%
filter(str_detect(ball, "ballA|ballB"),
str_detect(object, "ball")) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
arrange(clip, time) %>%
group_by(clip) %>%
summarize(A = any(!is.na(AE)),
B = any(!is.na(BE)),
A = ifelse(any(!is.na(AB)) &
any(!is.na(BE)) &
max(AB, na.rm = T) < min(BE, na.rm = T),
T,
A),
B = ifelse(any(!is.na(AB)) &
any(!is.na(AE)) &
max(AB, na.rm = T) < min(AE, na.rm = T),
T,
B)) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "transfer") %>%
arrange(clip, ball)
# collect features
df.features = df.features.initial %>%
filter(ball != "E") %>%
mutate(ball = factor(ball)) %>%
mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
speed = abs(velx) + abs(vely)) %>%
left_join(tmp.contactE) %>%
left_join(tmp.impactE) %>%
left_join(tmp.impactTotal) %>%
left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
left_join(tmp.movingE) %>%
left_join(tmp.exclusive) %>%
select(-c(appearance, velx, vely))
# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])set.seed(1)
df.plot = df.exp3.counterfactual.long %>%
left_join(df.exp3.model.counterfactual %>%
select(clip, ball, prediction) %>%
mutate(ball = as.factor(ball)),
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo,
by = "clip") %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
2,
colorindex),
colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
2,
colorindex),
colorindex = ifelse(model == "prediction", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "prediction"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
# width = 12,
# height = 6)# best fitting model
df.exp3.counterfactual.means %>%
left_join(df.exp3.model.counterfactual,
by = c("clip", "ball")) %>%
summarize(simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.86 | 19.84 |
# deterministic model
df.exp3.counterfactual.means %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_a, outcome_b) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "outcome") %>%
mutate(ball = factor(ball,
levels = c("outcome_a", "outcome_b"),
labels = c("B", "A"))),
by = c("clip", "ball")) %>%
mutate(outcome = outcome * 100) %>%
summarize(simulation_r = cor(rating_mean, outcome),
simulation_rmse = rmse(rating_mean, outcome)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.83 | 30.28 |
df.data = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# using whether-causation from the model
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether +
(1 + whether | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_w")
fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
(1 + whether + how | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_wh")
fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
(1 + whether + how + sufficient | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_whs")fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
elpd_diff se_diff
fit.brm_exp3_causal_whs 0.0 0.0
fit.brm_exp3_causal_wh -125.2 15.6
fit.brm_exp3_causal_w -426.2 31.1
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving)
# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving +
speed +
contact_e +
e_speed_diff +
e_direction_diff +
total_speed_diff +
total_direction_diff +
transfer +
e_moving +
exclusive + (1 | participant),
prior = priors,
data = df.regression.features %>%
select(-c(clip, ball, order, clipindex,
contains("outcome"))),
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant)
Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8.53 1.22 6.44 11.18 1.00 1163 1936
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 49.79 1.46 46.91 52.74 1.00 1129
moving 0.22 0.22 0.01 0.81 1.00 3399
speed 2.06 0.85 0.44 3.67 1.00 2319
contact_e 0.38 0.35 0.01 1.29 1.00 3201
e_speed_diff 0.12 0.12 0.00 0.46 1.00 4435
e_direction_diff 1.04 0.74 0.04 2.71 1.00 2202
total_speed_diff 2.21 0.95 0.41 4.17 1.00 2208
total_direction_diff 3.99 0.91 2.14 5.68 1.00 2718
transfer 15.59 0.79 14.06 17.14 1.00 3555
e_moving 0.18 0.17 0.01 0.62 1.00 3296
exclusive 4.39 0.73 2.97 5.81 1.00 3526
Tail_ESS
Intercept 1854
moving 1547
speed 1266
contact_e 2020
e_speed_diff 2268
e_direction_diff 1685
total_speed_diff 1393
total_direction_diff 1742
transfer 2992
e_moving 1738
exclusive 2538
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.21 0.46 32.33 34.10 1.00 5313 2845
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
func_fit_data = function(fit){
fit %>%
fitted(newdata = df.exp3.model.causal %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive),
.funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving),
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate)
}
df.exp3.causal.regression = df.exp3.causal.means %>%
mutate(w = func_fit_data(fit.brm_exp3_causal_w),
wh = func_fit_data(fit.brm_exp3_causal_wh),
whs = func_fit_data(fit.brm_exp3_causal_whs),
heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))prediction = "whs"
df.exp3.causal.regression %>%
summarize(simulation_r = cor(rating_mean, .[[prediction]]),
simulation_rmse = rmse(rating_mean, .[[prediction]])) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.87 | 13.07 |
df.fit = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# priors
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline =
brm(formula = rating ~ 1,
data = df.fit %>%
filter(participant == 1),
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_baseline"))
fit.brm_exp3_causal_individual_w =
brm(formula = rating ~ 1 + whether,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_w"))
fit.brm_exp3_causal_individual_wh =
brm(formula = rating ~ 1 + whether + how,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_wh"))
fit.brm_exp3_causal_individual_whs =
brm(formula = rating ~ 1 + whether + how + sufficient,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_whs"))
# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
group_by(participant) %>%
nest() %>%
ungroup() %>%
# fit model for each participant
mutate(fit_baseline = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_baseline,
newdata = .x)),
fit_w = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_w,
newdata = .x)),
fit_wh = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_wh,
newdata = .x)),
fit_whs = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_whs,
newdata = .x))) %>%
mutate(fit_baseline = map(.x = fit_baseline,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_w = map(.x = fit_w,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_wh = map(.x = fit_wh,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_whs = map(.x = fit_whs,
~ add_criterion(.x, criterion = c("loo", "waic"))),
r_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
model_comparison = pmap(.l = list(baseline = fit_baseline,
w = fit_w,
wh = fit_wh,
whs = fit_whs),
.f = ~ loo_compare(..1, ..2, ..3, ..4)),
best_model = map_chr(.x = model_comparison,
.f = ~ rownames(.) %>% .[1]),
best_model = factor(best_model,
levels = c("..1", "..2", "..3", "..4"),
labels = c("baseline", "w", "wh", "whs")))
save(list = c("df.exp3.causal.individual_fit"),
file = "data/exp3_causal_individual_fit.RData")load("data/exp3_causal_individual_fit.RData")
# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>%
count(best_model) %>%
print_table()| best_model | n |
|---|---|
| wh | 2 |
| whs | 39 |
df.exp3.causal.individual_fit %>%
select(participant, contains("r_"), contains("rmse_")) %>%
pivot_longer(cols = -participant,
names_to = c("measure", "model"),
names_sep = "_",
values_to = "fit") %>%
group_by(model, measure) %>%
summarize(quantiles = quantile(fit, probs = c(0.05, 0.5, 0.95)),
prob = c(0.05, 0.5, 0.95)) %>%
pivot_wider(names_from = prob,
values_from = quantiles) %>%
mutate(across(.fns = ~ round(., 2))) %>%
print_table()| model | measure | 0.05 | 0.5 | 0.95 |
|---|---|---|---|---|
| w | r | 0.18 | 0.43 | 0.57 |
| w | rmse | 29.33 | 36.03 | 43.62 |
| wh | r | 0.35 | 0.60 | 0.78 |
| wh | rmse | 23.82 | 32.19 | 40.60 |
| whs | r | 0.39 | 0.67 | 0.79 |
| whs | rmse | 22.69 | 31.17 | 38.12 |
set.seed(1)
df.cluster_whs = fit.brm_exp3_causal_whs %>%
ranef() %>%
.$participant %>%
as_tibble() %>%
select(contains("Estimate")) %>%
set_names(tolower(str_replace_all(names(.), "Estimate.", ""))) %>%
mutate(participant = 1:n()) %>%
relocate(participant)
df.cluster_whs = df.cluster_whs %>%
mutate(cluster = df.cluster_whs %>%
select(-participant) %>%
kmeans(centers = 2) %>%
.$cluster)
df.cluster_whs %>%
print_table()| participant | intercept | whether | how | sufficient | cluster |
|---|---|---|---|---|---|
| 1 | 8.30 | 0.44 | -5.36 | -0.84 | 1 |
| 2 | 0.03 | -20.89 | 32.67 | -20.83 | 2 |
| 3 | 2.27 | -8.27 | 26.27 | -9.42 | 2 |
| 4 | -0.66 | 7.98 | -2.96 | 10.74 | 1 |
| 5 | -0.75 | 2.47 | -5.92 | 0.20 | 1 |
| 6 | 5.09 | -15.71 | 22.42 | -18.09 | 2 |
| 7 | -0.17 | 1.71 | -3.13 | 0.22 | 1 |
| 8 | -3.55 | -2.49 | 4.88 | -0.21 | 2 |
| 9 | -4.15 | -5.31 | 9.38 | -2.57 | 2 |
| 10 | 2.10 | -9.99 | 3.25 | -10.96 | 2 |
| 11 | 2.10 | 7.99 | -10.81 | 5.96 | 1 |
| 12 | -2.00 | 15.83 | -25.89 | 18.27 | 1 |
| 13 | 3.03 | 4.98 | 0.72 | 2.29 | 1 |
| 14 | 3.17 | 4.02 | -11.54 | 1.56 | 1 |
| 15 | -5.07 | -13.56 | -3.49 | -11.90 | 2 |
| 16 | -0.89 | 23.98 | -17.31 | 25.30 | 1 |
| 17 | 0.35 | -4.04 | 2.84 | -3.79 | 2 |
| 18 | 3.05 | -0.15 | -10.15 | -2.03 | 1 |
| 19 | -2.30 | -7.53 | 4.17 | -7.76 | 2 |
| 20 | 8.44 | 20.15 | -11.67 | 17.52 | 1 |
| 21 | -1.69 | 1.54 | -12.89 | 2.37 | 1 |
| 22 | -5.38 | 1.33 | -3.82 | 2.56 | 1 |
| 23 | -0.53 | 1.24 | -7.19 | 3.29 | 1 |
| 24 | -0.43 | 1.91 | -0.12 | 2.55 | 1 |
| 25 | 2.88 | 12.33 | -12.26 | 9.65 | 1 |
| 26 | 0.46 | 12.55 | 1.40 | 11.74 | 1 |
| 27 | -4.46 | 16.00 | -20.07 | 21.43 | 1 |
| 28 | 2.21 | -8.24 | -1.57 | -10.35 | 2 |
| 29 | -0.78 | 7.61 | -4.96 | 7.93 | 1 |
| 30 | -3.78 | -21.34 | 18.53 | -19.98 | 2 |
| 31 | -6.59 | 7.50 | -12.71 | 12.17 | 1 |
| 32 | -1.70 | 2.50 | 3.37 | 6.52 | 1 |
| 33 | 3.29 | 13.81 | -10.56 | 13.26 | 1 |
| 34 | -4.20 | -12.18 | 9.83 | -10.61 | 2 |
| 35 | 0.41 | -16.26 | 0.37 | -18.28 | 2 |
| 36 | -1.55 | -27.69 | 47.62 | -28.75 | 2 |
| 37 | 10.64 | 13.67 | 0.43 | 8.39 | 1 |
| 38 | 1.35 | 9.88 | -7.13 | 9.38 | 1 |
| 39 | -5.71 | -10.20 | 7.58 | -10.99 | 2 |
| 40 | -6.18 | -7.24 | 8.60 | -6.13 | 2 |
| 41 | 1.97 | -2.12 | -1.61 | -0.67 | 1 |
# A tibble: 2 x 2
cluster n
<int> <int>
1 1 25
2 2 16
A 2-cluster solution yields a good result according to a number of different validation measures (see here for more details on the different measures).
tmp = df.cluster_whs %>%
select(-participant, -cluster)
fit = clValid(obj = as.matrix(tmp),
nClust = 2:10,
clMethods = c("kmeans"),
validation = c("internal", "stability"))Warning in clValid(obj = as.matrix(tmp), nClust = 2:10, clMethods =
c("kmeans"), : rownames for data not specified, using 1:nrow(data)
Clustering Methods:
kmeans
Cluster sizes:
2 3 4 5 6 7 8 9 10
Validation Measures:
2 3 4 5 6 7 8 9 10
kmeans APN 0.0991 0.1419 0.1438 0.1181 0.2214 0.1734 0.2094 0.1964 0.2269
AD 20.3268 15.8118 13.3966 11.5686 11.4495 10.0034 9.6513 9.2184 8.8807
ADM 7.3494 4.1130 4.3773 2.2086 3.7885 3.1164 3.5423 3.9327 3.9086
FOM 8.3657 7.0101 6.0226 5.4730 5.1633 4.9791 4.8082 4.8092 4.9224
Connectivity 4.4262 15.5893 20.1163 22.5329 25.0329 28.1413 34.4710 42.6353 44.9687
Dunn 0.1046 0.1626 0.1687 0.2509 0.2509 0.2550 0.1815 0.2524 0.2524
Silhouette 0.4691 0.3871 0.3868 0.3738 0.3522 0.3255 0.2988 0.2846 0.2774
Optimal Scores:
Score Method Clusters
APN 0.0991 kmeans 2
AD 8.8807 kmeans 10
ADM 2.2086 kmeans 5
FOM 4.8082 kmeans 8
Connectivity 4.4262 kmeans 2
Dunn 0.2550 kmeans 7
Silhouette 0.4691 kmeans 2
set.seed(1)
model_index = "whs"
# model predictions
model_prediction = list(fit.brm_exp3_causal_w,
fit.brm_exp3_causal_wh,
fit.brm_exp3_causal_whs) %>%
map_dfr(~ fitted(., df.exp3.causal.means %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")),
re_formula = NA) %>%
as_tibble()) %>%
mutate(ball = rep(c("A", "B"), n()/2),
clip = rep(rep(1:32, each = 2), 3),
model = rep(c("w", "wh", "whs"),
each = 64))
df.plot = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_both),
by = c("clip")) %>%
left_join(model_prediction %>%
filter(model == model_index) %>%
select(model = Estimate, clip, ball),
by = c("clip", "ball")) %>%
pivot_longer(cols = c(rating, model),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" &
outcome_both == 0, 1, colorindex),
colorindex = ifelse(model == "rating" &
outcome_both == 1, 2, colorindex),
colorindex = ifelse(model == "model", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "model"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
# width = 12,
# height = 6)check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(3, 7, 15, 16, 23)) %>%
group_by(clip, ball) %>%
summarize(mean = mean(rating),
low = smean.cl.boot(rating)[2],
high = smean.cl.boot(rating)[3]) %>%
left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
by = c("clip", "ball")) %>%
ungroup() %>%
pivot_longer(cols = c(mean, w, wh, whs),
names_to = "index",
values_to = "value") %>%
mutate_at(.vars = vars(low, high),
.funs = ~ ifelse(index == "mean", ., NA)) %>%
mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain",
"double prevention",
"joint causation",
"overdetermination",
"preemption")))
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
value = -10,
index = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
index = NA,
value = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = value,
group = index,
fill = index)) +
geom_bar(stat = "identity", color = "black",
position = position_dodge(0.9),
width = 0.9) +
geom_errorbar(mapping = aes(ymin = low, ymax = high),
width = 0,
size = 1,
position = position_dodge(0.9)) +
annotation_custom(grob = ball_a,
xmin = 0.5,
xmax = 1.5,
ymin = -25,
ymax = -7) +
annotation_custom(grob = ball_b,
xmin = 1.5,
xmax = 2.5,
ymin = -25,
ymax = -7) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_grid(~clip) +
scale_fill_grey(start = 1,
end = 0,
labels = c("mean rating ", expression(CSM[W] ~ " "),
expression(CSM[WH] ~ " "),
expression(CSM[WHS] ~ " "))) +
labs(y = "causal responsibility", fill = "") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(legend.text.align = 0,
text = element_text(size = 20),
panel.grid = element_blank(),
legend.position = "bottom",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
legend.text = element_text(size = 30),
strip.text = element_text(size = 30),
legend.spacing = unit(0.5, "cm"),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(t = 5, unit = "cm"),
plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)
# ggsave("../../figures/plots/exp3_selection_bars.pdf",
# width = 20,
# height = 10,
# device = cairo_pdf)func_scatterplot = function(model){
if(model == "heuristic"){
xlabel = "Heuristic"
}else{
xlabel = bquote(CSM[.(toupper(model))])
}
l.models = list(w = fit.brm_exp3_causal_w,
wh = fit.brm_exp3_causal_wh,
whs = fit.brm_exp3_causal_whs,
heuristic = fit.brm_exp3_causal_heuristic)
df.data = df.exp3.causal.means %>%
left_join(df.exp3.model.causal, by = c("clip", "ball")) %>%
left_join(df.regression.features, by = c("clip", "ball"))
df.plot = l.models[[model]] %>%
fitted(newdata = df.data,
re_formula = NA) %>%
as_tibble() %>%
clean_names() %>%
bind_cols(df.data)
p = ggplot(data = df.plot,
mapping = aes(x = estimate,
y = rating_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_smooth(aes(y = estimate,
ymin = q2_5,
ymax = q97_5),
stat = "identity",
color = "black") +
geom_linerange(size = 0.75,
mapping = aes(ymin = rating_low,
ymax = rating_high),
color = "gray50") +
geom_point(size = 2) +
scale_color_manual(values = c("black", "#e41a1c"),
guide = F) +
coord_fixed(xlim = c(0, 100),
ylim = c(0, 100)) +
labs(x = xlabel,
y = "responsibility judgment") +
annotate(geom = "text",
label = str_c(
"r = ", cor(df.plot$estimate, df.plot$rating_mean) %>%
round(2) %>%
as.character() %>%
str_sub(start = 2),
"\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>%
round(2)),
x = 5,
y = 95,
size = 6,
hjust = 0) +
scale_x_continuous(breaks = seq(0, 100, 25)) +
scale_y_continuous(breaks = seq(0, 100, 25)) +
theme(text = element_text(size = 20))
return(p)
}
# for creating and saving an individual scatter plot
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"
# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
# width = 5,
# height = 5)
list(func_scatterplot(model = "w"),
func_scatterplot(model = "wh"),
func_scatterplot(model = "whs"),
func_scatterplot(model = "heuristic")) %>%
wrap_plots(ncol = 2)check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(7, 23, 3, 15, 16)) %>%
mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))
df.plot = df.plot %>%
left_join(df.cluster_whs %>%
group_by(cluster) %>%
mutate(group = factor(cluster,
levels = 1:2,
labels = str_c("n = ", group_size(.)))) %>%
ungroup() %>%
select(participant, cluster, group),
by = "participant")
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
rating = -10,
group = NA,
participant = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
group = NA,
participant = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = participant,
shape = group)) +
geom_line(mapping = aes(linetype = "individual",
color = group),
size = 1,
alpha = 0.3) +
geom_point(mapping = aes(color = group),
size = 1,
alpha = 0.3) +
stat_summary(mapping = aes(group = group,
color = group,
linetype = "mean"),
fun = "mean",
geom = "line",
size = 1.5) +
stat_summary(data = df.plot %>% filter(cluster == 1),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
stat_summary(data = df.plot %>% filter(cluster == 2),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
annotation_custom(grob = ball_a,
xmin = 0.5,
xmax = 1.5,
ymin = -30,
ymax = -10) +
annotation_custom(grob = ball_b,
xmin = 1.5,
xmax = 2.5,
ymin = -30,
ymax = -10) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_wrap(~ clip,
ncol = 8) +
coord_cartesian(xlim = c(0.9, 2.1),
ylim = c(-5, 105),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25),
labels = seq(0, 100, 25)) +
scale_color_brewer(palette = "Set1",
guide = "none") +
labs(y = "causal responsibility rating",
linetype = "legend",
shape = "") +
theme_bw() +
guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
keywidth = unit(1.2, "cm")),
shape = guide_legend(override.aes = list(color = c("#377eb8",
"#e41a1c"),
shape = c(19, 19),
alpha = c(1, 1),
size = c(5, 5)))) +
theme(panel.grid = element_blank(),
text = element_text(size = 20),
legend.position = c(0.505, 0.23),
axis.title.x = element_blank(),
legend.text = element_text(size = 20),
legend.title = element_text(size = 24),
legend.background = element_rect(fill = "transparent"),
strip.text = element_text(size = 30),
axis.text.x = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.box = "horizontal",
legend.spacing.x = unit(0.1, "cm"),
panel.spacing.x = unit(1, "cm"),
plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))
print(p)
# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
# plot = p,
# width = 20,
# height = 8.5,
# device = cairo_pdf)df.exp3.model %>%
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(difference, whether, how, sufficient, robust) %>%
correlate() %>%
shave() %>%
fashion() %>%
print_table()| rowname | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|
| difference | |||||
| whether | .50 | ||||
| how | .79 | .27 | |||
| sufficient | .21 | .10 | .36 | ||
| robust | .43 | .93 | .24 | .20 |
fit.brm_exp3_causal_w %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_w") %>%
bind_rows(fit.brm_exp3_causal_wh %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_wh")) %>%
bind_rows(fit.brm_exp3_causal_whs %>%
tidy(effects = "fixed") %>%
mutate(model = "CSM_whs")) %>%
mutate(term = tolower(term)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
select(model, everything(), -effect, -component) %>%
print_table()| model | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|
| CSM_w | (intercept) | 31.68 | 1.61 | 28.41 | 34.87 |
| CSM_w | whether | 39.30 | 2.16 | 35.11 | 43.54 |
| CSM_wh | (intercept) | 10.29 | 1.56 | 7.30 | 13.41 |
| CSM_wh | whether | 28.40 | 2.68 | 23.18 | 33.56 |
| CSM_wh | how | 36.07 | 2.60 | 31.01 | 41.17 |
| CSM_whs | (intercept) | 10.43 | 1.59 | 7.30 | 13.50 |
| CSM_whs | whether | 27.93 | 2.60 | 22.80 | 33.09 |
| CSM_whs | how | 26.47 | 2.91 | 20.83 | 32.42 |
| CSM_whs | sufficient | 30.90 | 3.06 | 25.03 | 36.93 |
df.exp3.causal.individual_fit %>%
select_if(~ (is.numeric(.) | is.factor(.))) %>%
mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>%
select(participant, everything(), best_model) %>%
print_table()| participant | r_w | r_wh | r_whs | rmse_w | rmse_wh | rmse_whs | best_model |
|---|---|---|---|---|---|---|---|
| 1 | 0.29 | 0.38 | 0.48 | 30.00 | 29.03 | 27.74 | whs |
| 2 | 0.20 | 0.77 | 0.77 | 39.32 | 27.94 | 27.55 | whs |
| 3 | 0.37 | 0.78 | 0.79 | 38.90 | 28.00 | 26.95 | whs |
| 4 | 0.41 | 0.55 | 0.63 | 43.62 | 40.60 | 38.31 | whs |
| 5 | 0.45 | 0.53 | 0.54 | 39.23 | 37.32 | 36.81 | whs |
| 6 | 0.18 | 0.54 | 0.54 | 40.84 | 36.32 | 36.07 | whs |
| 7 | 0.49 | 0.61 | 0.64 | 32.96 | 29.93 | 29.15 | whs |
| 8 | 0.39 | 0.63 | 0.68 | 38.37 | 33.26 | 31.64 | whs |
| 9 | 0.40 | 0.72 | 0.76 | 35.86 | 28.34 | 26.39 | whs |
| 10 | 0.28 | 0.52 | 0.53 | 29.13 | 26.28 | 25.85 | whs |
| 11 | 0.51 | 0.56 | 0.60 | 33.66 | 32.19 | 31.17 | whs |
| 12 | 0.56 | 0.57 | 0.70 | 32.78 | 32.08 | 28.84 | whs |
| 13 | 0.56 | 0.71 | 0.74 | 29.83 | 25.40 | 24.21 | whs |
| 14 | 0.43 | 0.47 | 0.50 | 33.72 | 32.82 | 32.23 | whs |
| 15 | 0.19 | 0.35 | 0.37 | 39.77 | 38.29 | 37.94 | whs |
| 16 | 0.62 | 0.67 | 0.76 | 40.49 | 37.98 | 34.39 | whs |
| 17 | 0.42 | 0.66 | 0.71 | 28.57 | 23.82 | 22.42 | whs |
| 18 | 0.31 | 0.35 | 0.39 | 37.78 | 37.09 | 36.57 | whs |
| 19 | 0.37 | 0.60 | 0.62 | 33.41 | 29.23 | 28.62 | whs |
| 20 | 0.57 | 0.61 | 0.68 | 36.47 | 35.00 | 32.71 | whs |
| 21 | 0.43 | 0.50 | 0.57 | 30.95 | 29.61 | 28.32 | whs |
| 22 | 0.51 | 0.67 | 0.71 | 33.72 | 29.36 | 28.00 | whs |
| 23 | 0.33 | 0.45 | 0.52 | 38.76 | 37.00 | 35.57 | whs |
| 24 | 0.46 | 0.63 | 0.69 | 34.01 | 29.96 | 28.22 | whs |
| 25 | 0.61 | 0.66 | 0.70 | 30.60 | 28.91 | 27.55 | whs |
| 26 | 0.57 | 0.71 | 0.76 | 40.43 | 35.03 | 32.80 | whs |
| 27 | 0.46 | 0.52 | 0.66 | 43.12 | 41.49 | 38.12 | whs |
| 28 | 0.30 | 0.45 | 0.46 | 29.33 | 27.62 | 27.38 | whs |
| 29 | 0.48 | 0.60 | 0.65 | 38.93 | 36.03 | 34.48 | whs |
| 30 | 0.18 | 0.70 | 0.70 | 33.06 | 24.95 | 24.71 | whs |
| 31 | 0.48 | 0.59 | 0.69 | 37.54 | 34.72 | 31.76 | whs |
| 32 | 0.37 | 0.61 | 0.70 | 40.62 | 35.94 | 32.92 | whs |
| 33 | 0.50 | 0.56 | 0.64 | 38.64 | 36.79 | 34.85 | whs |
| 34 | 0.24 | 0.51 | 0.52 | 44.50 | 40.86 | 40.19 | whs |
| 35 | 0.15 | 0.32 | 0.30 | 34.56 | 33.36 | 33.37 | wh |
| 36 | 0.19 | 0.84 | 0.82 | 45.38 | 29.58 | 29.70 | wh |
| 37 | 0.55 | 0.63 | 0.67 | 36.04 | 33.36 | 32.11 | whs |
| 38 | 0.52 | 0.61 | 0.67 | 36.03 | 33.39 | 31.53 | whs |
| 39 | 0.47 | 0.77 | 0.77 | 29.61 | 21.75 | 21.58 | whs |
| 40 | 0.46 | 0.78 | 0.79 | 32.41 | 23.60 | 22.69 | whs |
| 41 | 0.32 | 0.50 | 0.58 | 32.90 | 30.38 | 28.69 | whs |
df.exp3.model.causal %>%
left_join(df.exp3.causal.regression,
by = c("clip", "ball")) %>%
filter(clip %in% c(7, 23, 16, 3, 15)) %>%
mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
.funs = ~ round(., 2)) %>%
select(clip, ball, difference, whether, how, sufficient, robust) %>%
mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ ifelse(. < 0.5,
str_c("xmark (", ., ")"),
str_c("cmark (", ., ")"))) %>%
arrange(clip) %>%
print_table()| clip | ball | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|---|
| 7 | A | cmark (1) | xmark (0.34) | cmark (1) | xmark (0) | xmark (0.25) |
| 7 | B | cmark (1) | cmark (1) | cmark (1) | cmark (0.67) | cmark (0.6) |
| 23 | A | xmark (0.05) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 23 | B | cmark (0.91) | cmark (0.79) | xmark (0) | xmark (0) | cmark (0.72) |
| 16 | A | cmark (1) | xmark (0.23) | cmark (1) | cmark (1) | xmark (0.35) |
| 16 | B | xmark (0) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 3 | A | cmark (1) | cmark (0.88) | cmark (1) | xmark (0.12) | cmark (0.76) |
| 3 | B | cmark (1) | cmark (0.89) | cmark (1) | xmark (0.11) | cmark (0.75) |
| 15 | A | cmark (1) | xmark (0.01) | cmark (1) | cmark (0.99) | xmark (0.1) |
| 15 | B | cmark (1) | xmark (0.01) | cmark (1) | cmark (1) | xmark (0.1) |
df.exp3.causal.regression %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo %>%
select(-clipindex),
by = c("clip")) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ . * 100) %>%
mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
select(clip, ball, contains("outcome"), difference, whether, how, sufficient,
robust, w, wh, whs, heuristic, rating = rating_mean) %>%
# print_table(format = "latex", digits = 0)
print_table(digits = 0)| clip | ball | outcome_both | outcome_a | outcome_b | outcome_none | difference | whether | how | sufficient | robust | w | wh | whs | heuristic | rating |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | A | 0 | 0 | 0 | 0 | 100 | 40 | 100 | 23 | 36 | 47 | 58 | 55 | 57 | 42 |
| 1 | B | 0 | 0 | 0 | 0 | 100 | 15 | 100 | 16 | 9 | 38 | 51 | 46 | 54 | 37 |
| 2 | A | 0 | 0 | 0 | 0 | 57 | 12 | 0 | 0 | 10 | 36 | 14 | 14 | 25 | 21 |
| 2 | B | 0 | 0 | 0 | 0 | 18 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 24 | 19 |
| 3 | A | 1 | 0 | 0 | 0 | 100 | 88 | 100 | 12 | 76 | 66 | 71 | 65 | 72 | 76 |
| 3 | B | 1 | 0 | 0 | 0 | 100 | 89 | 100 | 11 | 75 | 67 | 72 | 65 | 72 | 75 |
| 4 | A | 1 | 0 | 0 | 0 | 100 | 78 | 100 | 4 | 78 | 62 | 69 | 60 | 58 | 63 |
| 4 | B | 1 | 0 | 0 | 0 | 100 | 95 | 100 | 15 | 57 | 69 | 73 | 68 | 54 | 78 |
| 5 | A | 0 | 0 | 1 | 0 | 100 | 90 | 100 | 0 | 47 | 67 | 72 | 62 | 47 | 71 |
| 5 | B | 0 | 0 | 1 | 0 | 100 | 0 | 100 | 0 | 0 | 32 | 46 | 37 | 68 | 22 |
| 6 | A | 0 | 0 | 1 | 0 | 100 | 59 | 100 | 16 | 35 | 55 | 63 | 59 | 53 | 73 |
| 6 | B | 0 | 0 | 1 | 0 | 100 | 18 | 100 | 6 | 14 | 39 | 52 | 44 | 53 | 22 |
| 7 | A | 1 | 0 | 1 | 0 | 100 | 34 | 100 | 0 | 25 | 45 | 56 | 46 | 70 | 59 |
| 7 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 67 | 60 | 71 | 75 | 86 | 64 | 79 |
| 8 | A | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 25 | 7 |
| 8 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 100 | 100 | 71 | 75 | 96 | 84 | 92 |
| 9 | A | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 8 |
| 9 | B | 0 | 1 | 0 | 0 | 100 | 100 | 100 | 0 | 100 | 71 | 75 | 65 | 78 | 90 |
| 10 | A | 0 | 1 | 0 | 0 | 77 | 18 | 0 | 0 | 22 | 39 | 15 | 16 | 15 | 23 |
| 10 | B | 0 | 1 | 0 | 0 | 98 | 79 | 0 | 0 | 63 | 63 | 33 | 33 | 21 | 55 |
| 11 | A | 1 | 1 | 0 | 0 | 100 | 70 | 100 | 77 | 68 | 59 | 66 | 80 | 71 | 93 |
| 11 | B | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 4 |
| 12 | A | 1 | 1 | 0 | 0 | 100 | 82 | 100 | 74 | 83 | 64 | 70 | 83 | 57 | 77 |
| 12 | B | 1 | 1 | 0 | 0 | 100 | 0 | 100 | 12 | 24 | 32 | 46 | 41 | 53 | 37 |
| 13 | A | 0 | 1 | 1 | 0 | 67 | 34 | 0 | 0 | 35 | 45 | 20 | 20 | 14 | 8 |
| 13 | B | 0 | 1 | 1 | 0 | 70 | 35 | 0 | 0 | 35 | 46 | 20 | 20 | 21 | 64 |
| 14 | A | 0 | 1 | 1 | 0 | 97 | 91 | 0 | 0 | 59 | 67 | 36 | 36 | 21 | 22 |
| 14 | B | 0 | 1 | 1 | 0 | 91 | 77 | 0 | 0 | 51 | 62 | 32 | 32 | 20 | 18 |
| 15 | A | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 99 | 10 | 32 | 47 | 68 | 71 | 76 |
| 15 | B | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 100 | 10 | 32 | 47 | 68 | 71 | 76 |
| 16 | A | 1 | 1 | 1 | 0 | 100 | 23 | 100 | 100 | 35 | 41 | 53 | 74 | 80 | 92 |
| 16 | B | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 4 |
| 17 | A | 0 | 0 | 0 | 1 | 100 | 19 | 100 | 37 | 18 | 39 | 52 | 54 | 66 | 69 |
| 17 | B | 0 | 0 | 0 | 1 | 100 | 0 | 100 | 36 | 17 | 32 | 46 | 48 | 65 | 46 |
| 18 | A | 0 | 0 | 0 | 1 | 100 | 11 | 100 | 40 | 17 | 36 | 49 | 52 | 55 | 63 |
| 18 | B | 0 | 0 | 0 | 1 | 100 | 7 | 100 | 37 | 9 | 34 | 48 | 50 | 56 | 66 |
| 19 | A | 1 | 0 | 0 | 1 | 100 | 74 | 100 | 7 | 65 | 61 | 67 | 60 | 55 | 53 |
| 19 | B | 1 | 0 | 0 | 1 | 100 | 72 | 100 | 7 | 65 | 60 | 67 | 59 | 55 | 49 |
| 20 | A | 1 | 0 | 0 | 1 | 100 | 92 | 100 | 8 | 72 | 68 | 72 | 65 | 57 | 41 |
| 20 | B | 1 | 0 | 0 | 1 | 100 | 88 | 100 | 4 | 53 | 66 | 71 | 63 | 56 | 71 |
| 21 | A | 0 | 0 | 1 | 1 | 100 | 47 | 100 | 40 | 45 | 50 | 60 | 63 | 58 | 80 |
| 21 | B | 0 | 0 | 1 | 1 | 100 | 9 | 100 | 21 | 10 | 35 | 49 | 46 | 59 | 18 |
| 22 | A | 0 | 0 | 1 | 1 | 100 | 100 | 100 | 89 | 83 | 71 | 75 | 92 | 48 | 60 |
| 22 | B | 0 | 0 | 1 | 1 | 100 | 8 | 100 | 0 | 15 | 35 | 49 | 39 | 53 | 42 |
| 23 | A | 1 | 0 | 1 | 1 | 5 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 15 | 3 |
| 23 | B | 1 | 0 | 1 | 1 | 91 | 79 | 0 | 0 | 72 | 63 | 33 | 33 | 22 | 39 |
| 24 | A | 1 | 0 | 1 | 1 | 100 | 66 | 100 | 4 | 63 | 58 | 65 | 57 | 57 | 44 |
| 24 | B | 1 | 0 | 1 | 1 | 100 | 94 | 100 | 22 | 79 | 69 | 73 | 70 | 54 | 73 |
| 25 | A | 0 | 1 | 0 | 1 | 100 | 25 | 100 | 21 | 26 | 42 | 53 | 50 | 69 | 43 |
| 25 | B | 0 | 1 | 0 | 1 | 100 | 74 | 100 | 54 | 65 | 61 | 67 | 74 | 56 | 73 |
| 26 | A | 0 | 1 | 0 | 1 | 100 | 6 | 100 | 3 | 9 | 34 | 48 | 40 | 60 | 39 |
| 26 | B | 0 | 1 | 0 | 1 | 100 | 87 | 100 | 35 | 54 | 66 | 71 | 72 | 46 | 69 |
| 27 | A | 1 | 1 | 0 | 1 | 100 | 97 | 100 | 52 | 97 | 70 | 74 | 80 | 67 | 80 |
| 27 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 17 | 6 |
| 28 | A | 1 | 1 | 0 | 1 | 100 | 90 | 100 | 22 | 80 | 67 | 72 | 69 | 79 | 89 |
| 28 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 12 | 5 |
| 29 | A | 0 | 1 | 1 | 1 | 100 | 58 | 100 | 24 | 44 | 55 | 63 | 61 | 66 | 47 |
| 29 | B | 0 | 1 | 1 | 1 | 100 | 63 | 100 | 24 | 38 | 57 | 64 | 62 | 54 | 67 |
| 30 | A | 0 | 1 | 1 | 1 | 100 | 57 | 100 | 29 | 49 | 54 | 63 | 62 | 63 | 58 |
| 30 | B | 0 | 1 | 1 | 1 | 100 | 46 | 100 | 24 | 39 | 50 | 60 | 57 | 63 | 56 |
| 31 | A | 1 | 1 | 1 | 1 | 100 | 2 | 100 | 4 | 3 | 32 | 47 | 39 | 63 | 44 |
| 31 | B | 1 | 1 | 1 | 1 | 100 | 4 | 100 | 4 | 4 | 33 | 47 | 39 | 53 | 46 |
| 32 | A | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 5 |
| 32 | B | 1 | 1 | 1 | 1 | 100 | 75 | 100 | 66 | 73 | 61 | 68 | 78 | 65 | 71 |
fit.brm_exp3_causal_heuristic %>%
tidy(effects = "fixed") %>%
mutate(term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = conf.low,
`upper 95% HDI` = conf.high) %>%
print_table()Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
| effect | component | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|---|
| fixed | cond | (intercept) | 49.79 | 1.46 | 46.91 | 52.74 |
| fixed | cond | moving | 0.22 | 0.22 | 0.01 | 0.81 |
| fixed | cond | speed | 2.06 | 0.85 | 0.44 | 3.67 |
| fixed | cond | contact_e | 0.38 | 0.35 | 0.01 | 1.29 |
| fixed | cond | e_speed_diff | 0.12 | 0.12 | 0.00 | 0.46 |
| fixed | cond | e_direction_diff | 1.04 | 0.74 | 0.04 | 2.71 |
| fixed | cond | total_speed_diff | 2.21 | 0.95 | 0.41 | 4.17 |
| fixed | cond | total_direction_diff | 3.99 | 0.91 | 2.14 | 5.68 |
| fixed | cond | transfer | 15.59 | 0.79 | 14.06 | 17.14 |
| fixed | cond | e_moving | 0.18 | 0.17 | 0.01 | 0.62 |
| fixed | cond | exclusive | 4.39 | 0.73 | 2.97 | 5.81 |
We have run a pilot eye-tracking experiment with 15 participants who saw the same clips as the ones in Experiment 3. Instead of asking participants to judge “To what extent were A and B responsible for E (not) going through the gate?”, we asked them to judge to what extent they agreed with the following statement: “Ball A/B caused ball E to go through the gate.” or “Ball A/B prevented ball E from going through the gate.” depending on the outcome.
Participants’ agreement judgments with the causal statements in this pilot eye-tracking experiment were strikingly similar to participants’ responsibility judgments in the online experiment.
df.exp3.eye_tracking_pilot =
read_csv("../../data/experiment3_eye_tracking_pilot.csv") %>%
rename(clip = trial) %>%
mutate(ball = str_to_upper(ball)) %>%
group_by(clip, ball, outcome) %>%
summarize(causality_mean = mean(causation),
causality_low = smean.cl.boot(causation)[2],
causality_high = smean.cl.boot(causation)[3]) %>%
ungroup()
df.plot = df.exp3.eye_tracking_pilot %>%
# select(clip, ball, outcome,) %>%
left_join(df.exp3.causal.means %>%
select(clip,
ball,
responsibility_mean = rating_mean,
responsibility_low = rating_low,
responsibility_high = rating_high),
by = c("clip", "ball"))
# highlight the double prevention clip (#23)
df.plot = df.plot %>%
mutate(color = ifelse(clip == 23, "1", "0"))
ggplot(data = df.plot,
mapping = aes(x = causality_mean,
y = responsibility_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_linerange(mapping = aes(ymin = responsibility_low,
ymax = responsibility_high),
alpha = 0.1) +
geom_linerange(mapping = aes(xmin = causality_low,
xmax = causality_high),
alpha = 0.1) +
geom_point(size = 2,
mapping = aes(color = color),
show.legend = F) +
geom_smooth(method = "lm",
color = "black") +
annotate(geom = "text",
x = c(0, 0),
y = c(100, 90),
label = c(str_c("r = ",
cor(df.plot$causality_mean,
df.plot$responsibility_mean) %>%
round(2)),
str_c("RMSE = ", rmse(df.plot$causality_mean,
df.plot$responsibility_mean) %>%
round(2))),
hjust = 0,
size = 8) +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_color_manual(values = c("black", "red")) +
labs(x = "causal judgment",
y = "responsibility judgment") +
theme(text = element_text(size = 24))
# ggsave("../../figures/plots/aux_question_framing_comparison.pdf",
# width = 6,
# height = 6)df.plot = fit.brm_exp3_causal_whs %>%
fitted() %>%
as_tibble() %>%
select(prediction = Estimate) %>%
bind_cols(df.exp3.causal.long) %>%
relocate(prediction, .after = last_col())
df.text = df.plot %>%
group_by(participant) %>%
summarize(r = round(cor(prediction, rating), 2)) %>%
ungroup() %>%
mutate(r = str_c("r = ", r),
prediction = 1,
rating = 110)
ggplot(data = df.plot,
mapping = aes(x = prediction,
y = rating)) +
geom_abline(intercept = 0,
slope = 1,
color = "blue",
alpha = 0.5) +
geom_point(alpha = 0.3,
size = 1) +
geom_text(data = df.text,
mapping = aes(label = r),
hjust = 0,
color = "red",
size = 4) +
facet_wrap(facets = vars(participant),
ncol = 7) +
labs(x = "model prediction",
y = "causal responsibility rating") +
coord_cartesian(xlim = c(0, 100),
ylim = c(0, 100),
expand = F,
clip = "off") +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
expand = expansion(mult = c(0, 0))) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_blank(),
text = element_text(size = 16),
panel.spacing.x = unit(0.75, "cm"),
panel.spacing.y = unit(1, "cm"),
plot.margin = margin(t = 0.7,
r = 0.8,
b = 0.2,
l = 0.2,
unit = "cm"))
# ggsave("../../figures/plots/exp3_individual_scatter_plots.pdf",
# width = 12,
# height = 8)df.plot = df.exp3.causal.individual_fit %>%
select(participant, contains("r_"), contains("rmse_")) %>%
pivot_longer(cols = -participant,
names_to = c("measure", "model"),
names_sep = "_",
values_to = "fit")
ggplot(data = df.plot %>%
filter(measure == "r"),
mapping = aes(x = fit,
fill = model)) +
geom_density(alpha = 0.5) +
labs(x = "correlation value") +
scale_x_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_fill_brewer(palette = "Set1") +
# ggplot2::theme_classic() +
theme(legend.position = c(0.3, 0.95),
legend.direction = "horizontal",
text = element_text(size = 20))
# ggsave("../../figures/plots/exp3_individual_densities.pdf",
# width = 8,
# height = 6)
### 3D scatter plot of participant clusters
plot_ly(x = df.cluster_whs$whether,
y = df.cluster_whs$how,
z = df.cluster_whs$sufficient,
type = "scatter3d",
mode = "markers",
color = as.factor(df.cluster_whs$cluster),
colors = c("#e41a1c", "#377eb8")) %>%
layout(showlegend = FALSE,
title = "Participant clusters",
scene = list(
xaxis = list(title = "whether"),
yaxis = list(title = "how"),
zaxis = list(title = "sufficient")))set.seed(1)
df.plot = df.exp3.causal.long %>%
left_join(df.cluster_whs %>%
select(participant, cluster) %>%
mutate(participant = as.numeric(participant),
cluster = as.factor(cluster)) %>%
group_by(cluster) %>%
mutate(label = factor(cluster,
levels = 1:2,
labels = str_c("n = ", group_size(.)))) %>%
ungroup(),
by = "participant")
ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = label,
fill = label)) +
geom_point(shape = 21,
alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
jitter.height = 0,
dodge.width = 0.75)) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
mapping = aes(color = label),
size = 1,
position = position_dodge(width = 0.75)) +
stat_summary(fun = "mean",
geom = "point",
size = 4,
shape = 21,
position = position_dodge(width = 0.75)) +
facet_wrap(facets = vars(clip), ncol = 8) +
labs(y = "causal responsibility rating",
fill = "cluster",
color = "cluster") +
scale_fill_brewer(palette = "Set1") +
scale_color_brewer(palette = "Set1") +
theme(text = element_text(size = 24),
legend.position = "bottom")
# ggsave("../../figures/plots/exp3_clusters_points.pdf",
# width = 12,
# height = 8)library("ggtern")
df.plot = df.exp3.causal.individual_fit %>%
select(participant, fit_whs) %>%
mutate(estimates = map(fit_whs, ~ fixef(.) %>%
as_tibble(rownames = "term") %>%
clean_names())) %>%
select(participant, estimates) %>%
unnest(estimates) %>%
filter(!str_detect(term, "Intercept")) %>%
select(participant, term, estimate) %>%
pivot_wider(names_from = term,
values_from = estimate) %>%
mutate_at(vars(whether, how, sufficient),
.funs = list(norm = ~ . / (how + whether + sufficient))) %>%
mutate(color = 0,
color = ifelse(how_norm == max(how_norm), 1, color),
color = ifelse(whether_norm == max(whether_norm), 2, color),
color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
color = factor(color))
ggplot(data = df.plot,
mapping = aes(x = how,
y = sufficient,
z = whether,
color = color)) +
geom_point(alpha = 0.7,
size = 2,
show.legend = F) +
scale_color_manual(values = c("black", "red", "green", "blue")) +
coord_tern() +
theme_showarrows() +
theme(text = element_text(size = 20),
tern.axis.title.T = element_blank(),
tern.axis.title.L = element_blank(),
tern.axis.title.R = element_blank())
# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
# width = 5,
# height = 5)R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.0 purrr_0.3.4
[5] readr_1.3.1 tidyr_1.1.0 tibble_3.0.1 tidyverse_1.3.0
[9] patchwork_1.0.0 clValid_0.6-7 cluster_2.1.0 egg_0.4.5
[13] gridExtra_2.3 plotly_4.9.2.1 png_0.1-7 tidybayes_2.1.1
[17] Hmisc_4.4-0 ggplot2_3.3.1 Formula_1.2-3 survival_3.1-8
[21] lattice_0.20-38 janitor_2.0.1 jsonlite_1.6.1 xtable_1.8-4
[25] corrr_0.4.2 brms_2.13.0 Rcpp_1.0.4.6 broom.mixed_0.2.6
[29] lme4_1.1-23 Matrix_1.2-18 kableExtra_1.1.0 knitr_1.28
loaded via a namespace (and not attached):
[1] readxl_1.3.1 backports_1.1.7 plyr_1.8.6
[4] igraph_1.2.5 lazyeval_0.2.2 svUnit_1.0.3
[7] TMB_1.7.16 splines_3.6.3 crosstalk_1.1.0.1
[10] rstantools_2.1.0 inline_0.3.15 digest_0.6.25
[13] htmltools_0.4.0 rsconnect_0.8.16 fansi_0.4.1
[16] magrittr_1.5 checkmate_2.0.0 modelr_0.1.8
[19] RcppParallel_5.0.1 matrixStats_0.56.0 xts_0.12-0
[22] prettyunits_1.1.1 jpeg_0.1-8.1 colorspace_1.4-1
[25] blob_1.2.1 rvest_0.3.5 ggdist_2.1.1
[28] haven_2.3.1 xfun_0.14 callr_3.4.3
[31] crayon_1.3.4 zoo_1.8-8 glue_1.4.1
[34] gtable_0.3.0 emmeans_1.4.7 webshot_0.5.2
[37] pkgbuild_1.0.8 rstan_2.19.3 abind_1.4-5
[40] scales_1.1.1 mvtnorm_1.1-1 DBI_1.1.0
[43] miniUI_0.1.1.1 viridisLite_0.3.0 htmlTable_1.13.3
[46] HDInterval_0.2.2 foreign_0.8-75 stats4_3.6.3
[49] StanHeaders_2.21.0-5 DT_0.13 htmlwidgets_1.5.1
[52] httr_1.4.1 threejs_0.3.3 arrayhelpers_1.1-0
[55] RColorBrewer_1.1-2 acepack_1.4.1 ellipsis_0.3.1
[58] farver_2.0.3 pkgconfig_2.0.3 loo_2.2.0
[61] dbplyr_1.4.4 nnet_7.3-12 utf8_1.1.4
[64] labeling_0.3 tidyselect_1.1.0 rlang_0.4.6
[67] reshape2_1.4.4 later_1.1.0.1 cellranger_1.1.0
[70] munsell_0.5.0 tools_3.6.3 cli_2.0.2
[73] generics_0.0.2 broom_0.7.0.9000 ggridges_0.5.2
[76] evaluate_0.14 fastmap_1.0.1 yaml_2.2.1
[79] fs_1.4.1 processx_3.4.2 nlme_3.1-144
[82] mime_0.9 xml2_1.3.2 compiler_3.6.3
[85] bayesplot_1.7.2 shinythemes_1.1.2 rstudioapi_0.11
[88] reprex_0.3.0 statmod_1.4.34 stringi_1.4.6
[91] highr_0.8 ps_1.3.3 Brobdingnag_1.2-6
[94] nloptr_1.2.2.1 markdown_1.1 shinyjs_1.1
[97] vctrs_0.3.1 pillar_1.4.4 lifecycle_0.2.0
[100] bridgesampling_1.0-0 estimability_1.3 data.table_1.12.8
[103] httpuv_1.5.4 R6_2.4.1 latticeExtra_0.6-29
[106] bookdown_0.19 promises_1.1.1 boot_1.3-24
[109] colourpicker_1.0 MASS_7.3-51.5 gtools_3.8.2
[112] assertthat_0.2.1 withr_2.2.0 shinystan_2.5.0
[115] mgcv_1.8-31 parallel_3.6.3 hms_0.5.3
[118] rpart_4.1-15 class_7.3-15 coda_0.19-3
[121] minqa_1.2.4 rmarkdown_2.2 snakecase_0.11.0
[124] shiny_1.4.0.2 lubridate_1.7.9 base64enc_0.1-3
[127] dygraphs_1.1.1.6